One-fluid description of turbulently flowing suspensions 
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We suggested a one-fluid model of a turbulent dilute suspension which accounts for the "two-way" 
fluid-particle interactions by fc-dependent effective density of suspension and additional damping 
term in the Navier-Stokes equation. We presented analytical description of the particle modifica- 
tion of turbulence including scale invariant suppression of a small k part of turbulent spectrum 
(independent of the particle response time) and possible enhancemenent of large k region [up to the 
factor (1 + <f>) 2 / 3 ]. Our results are in agreement with qualitative picture of isotropic homogeneous 
turbulence of dilute suspensions previously observed in laboratory and numerical experiments. 

PACS numbers: 47.57.Kf, 47.27.Gs, 47.10. +g 



Introduction. The interaction of solid particles or liq- 
uid droplets with turbulence in gases controls the perfor- 
mance of various engineering devices like the combustion 
of pulverized coal and liquid sprays and cyclone separa- 
tions. This interaction plays an important role in many 
areas of environmental science and physics of the atmo- 
sphere. Dust storms, rain triggering, dusting and spray- 
ing for agricultural or forestry purposes, preparation and 
processing of aerosols are typical examples. For a review 
of turbulent flows with particles and droplets see, e.g. jj]. 

In dilute suspensions with small volume fractions of 
particles, C p , the particle-particle interactions are negli- 
gible. Nevertheless, for p p /pi ~> 1 (the ratio of the solid 
particle material and the gas densities), the mass loading 
4> = CpPp/pf may exceed unity and kinetic energies of the 
particles and the carrier gas may be compatible; hence 
the "two-way coupling": effect of fluid on particles and 
vice versa must be accounted for. Current understanding 
of turbulence in dilute suspensions is still at its infancy 
due to the highly nonlinear nature of physically relevant 
interactions and a wide spectrum of governing parame- 
ters (the particle size a vs. L and 77, the outer and inner 
scales of turbulence, the particle response time r p vs. j L 
and 7,7, turnover frequencies of L- and 77- scale eddies). 

Analytical study of the problem is mainly based upon a 
two-fluid model wherein both the carrying fluid and par- 
ticle phases are treated as interpenetrating continua |l], 
||, H . This model deals with non-interacting solid spher- 
ical particles which are small enough such that: (i) one 
can neglect the effect of preferential concentration and 
assume homogeneity of the particle space distribution; 
(ii) the Stokes viscous drag law for particle acceleration, 
du p /dt — [v,f — m p ]/t p , is valid (ttf is the fluid veloc- 
ity). Unfortunately, statistical description of two-fluid 
turbulence by closure procedures requires a set of addi- 
tional questionable simplifications due to the lack of un- 
derstanding of the relevant physics of the particle-fluid 
interactions. This made closures of the two-fluid model 
highly qualitative at best ||, ||, |j| . 

We think that the basic physics of the problem may 
be described by a simpler one-fluid model of turbulent 



dilute suspensions, which requires standard closures of 
one-phase turbulence. The present Letter suggests such 
a model and, as a first step, uses a properly modified 
simple closure, based on the Kolmogorov-Richardson cas- 
cade picture of turbulence. The resulting non-linear dif- 
ferential equations for the energy budget were solved 
analytically with required accuracy. This provides an 
economical and internally consistent analytical descrip- 
tion of the turbulence modification by particles includ- 
ing the dependencies of suppression or enhancement of 
the turbulent activity on the three governing parameters: 
(T P 7 t ), 4> an d scale of eddies. These effects were previ- 
ously observed in numerous experimental and numerical 
studies |Q 0, |L |l H ||, but they still await analyti- 
cal description. Our analytical findings are in a qualita- 
tive agreements with the observations. We believe that 
the one-fluid model, together with more elaborated clo- 
sures of one-phase turbulence M , offers an insightful and 
constructive look at basic physics of even more complex 
particle-laden turbulent flows. 

1. One-fluid model of turbulent suspensions reads: 

r d 1 

Peff(fc) [g- + 7 P (fc)J u(t, k) + pk 2 u(t, k) (1) 
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Here u(t, k) is the incompressible velocity field of the car- 
rier fluid in the k representation and p is the dynamical 
viscosity. This equation differs from the Navier-Stokes 
(NS) equation in the three aspects: 

a. The carrier fluid density pf is replaced by p e s(k), a 
fc-dependent effective density of the suspension for tur- 
bulent fluctuations of scale 1/k [referred to as fc-eddies]; 

b. The fluid-particle friction is described by a damping 
frequency 7 P (fc). 

c. The interaction amplitude of the NS equation, 
7fc£fc 2 = Pt [P af} (k) W + P a ^{k)kP]6{k + fei + fc 2 )/2, 
P a P(k) = S a0 - k a kP/k 2 is replaced by 
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An underlying physics of these aspects is quite simple: 

a. For small k the turnover frequency j(k) of fc-eddies 
is small in the sense 7(fc)r p <C 1. Therefore in this region 
of k the particle velocity is very close to that of the carrier 
fluid and we can describe the suspension as one fluid with 
effective density p c s(k) which is close to the density of 
suspension, p s = pf + C p p p . In contrary, for large k, 
when 7(fc)r p 3> 1, the particles cannot follow these very 
fast motions and may be considered at rest. Thus the 
particles do not contribute to the effective density and 
Pcs(k) — > pt- In general case p c g(k) may be given by 



Pes(k) = pi [l + 0/com (fc)] ) = C P Pp/Pt 



(3) 



Here a statistical ensemble of all particles, partially in- 
volved in the motion of /c-eddies, is replaced by two sub- 
ensembles of "fully comoving" fraction f com (k) of parti- 
cles and "fully resting" fraction / res t(&) = 1 — /com (A) of 
particles, which does not contribute to p c e(k). 

b. The particles at rest cause fluid-particle friction. 
According to Newton's third law, the damping frequency 
of suspension 7 P (A:) may be related to the particle re- 
sponse time, r p , via the ratio of total mass of particles 
M p to the total effective mass of the suspension M e ff(k): 



7p(*) 



CpPp/rest(fc) _ 4>Pl /rcst(fc) 



, M off (k) 



, Pcs(k) 



(4) 



As we mentioned, the portions / com (fc) and f res t(k) de- 
pend on r p 7(fc). Clearly, for r p j(k) <C 1, f Te st{k) has 
the same smallness: f TCS t(k) ~ T p j(k). In the opposite 
case when l/r p 7(fc) is small, / com (k) has corresponding 
smallness: f CO m{k) ~ l/r p 7(fc). For r p 7(fc) = 1 one ex- 
pects f Ies t(k) ~ fcom(k) — \. As a simple model of such 
a function we adopt 

/ rest (fc) - l-/cam(*)=T P 7(fc)/[l + Tp7(fc)] J (5) 

which also follows from more elaborated analysis 0. 
With Eq. (||) we rewrite Eqs. (||) and (||) as follows: 

Peff(fc) = p f {l + 0/[l + r p7 (fc)]}, (6) 
7p (fc) = 7 (fc)/[l + + T- p7 (fc)] . (7) 

c. It may be shown that interaction amplitude ^kkk 
is Galilean invariant and conserves the total kinetic en- 
ergy of suspension (i.e. of the carrier fluid and fully co- 
moving part of the particles) if one neglects the damping 
terms pk 2 and 7 P (&) in Eq. (|l]). Notice, that the detailed 
form of fc is not important for this discussion, only 
the conservation of the energy is essential. 

The basic equations of our model, (|l|), (0), (||) and 
(Q) are self-consistent equations for turbulent velocity of 
the carrier fluid in which the coefficients depend on the 
eddy-turnover frequency "f(k) which, in its turn, depends 
on a stochastic solution of the same equations. 



2. Budget of the kinetic energy in suspensions. 

The definition of one-dimensional density of kinetic en- 
ergy in a single phase flow E(t, k) reads 



E(t,k) = p{k 2 F 2 (t,k)/2TT 



(8) 



Here F 2 (t, k) is the simultaneous pair velocity correla- 
tion for isotropic turbulence. Corresponding definition 
of one-dimensional density of kinetic energy for suspen- 
sions, £{k), has to account for the fc-dependent density: 



£(t,k) = p cS (k)k 2 F 2 (t,k)/2ir 



(9) 



Multiplying Eq. (|j) by u(t, k') and averaging, one gets 
the equation for the energy budget in the inertial interval: 

£(t,k)/2dt + lp (k)£(t,k) +de(k)/dk = . (10) 

Here e(k) is one- dimensional energy flux in the fc-space 
and we neglected the energy input in the outer scale L 
and the viscous damping, pk 2 , which becomes essential 
near the viscous microscale rj. 

In the Richardson-Kolmogorov picture of turbulence 
the only relevant parameters in the inertial interval are k, 
Pf and e(fc). Similarly, in our model ([!]) these parameters 
are k, p e s(k) and e(k). Other functions in the problem 
may be related to them by the dimensional reasoning: 

£(k) = Ci [e(k) 2 Pcff (k)] 1/3 k- 5 / 3 , (11) 

(12) 



7 (fc) = C 2 [e(k)/ PcS (k)] 1/3 k 2 ^. 



where C± ~ C 2 ~ 1 are dimensionless constants for par- 
ticle free case. Hereafter we omit the explicit reference 
to the time dependence. Substituting Eqs. (||), (@), (0) 
and (]l|) into Eq. (|o|) yields in the stationary case 



de(Jfe) e(fc) 



C» 



dk 



k l + + 7(fc)r p 



, C = Ci C 2 . (13) 



To find the iterative solution of Eq. (|i3| ) we denote as 
e n (k), Jn(k) and p e s, n (k) corresponding functions at nth 
iteration step; take for "zeroth step" their values at 
Pouter = l/L: e (k) = e(l/L) = e L , 70(A) = j{l/L) = 
7z,j Poff,o(^) = PeffiX/L) = p L ; define dimensionless func- 
tions of k = kL: s„(k) = e L f n (n), y n (k) = 7x,fl , n(«0, 
PeS,n(k) = Pfr n (n) and iterate the equations 



/n(«0 



ffn(«) 



exp 



-Ada; 



! x[l + Sg n -x(x)] 



r n -i(K) 

r n («) = 1 + (/>/[! + t p j L g n {K)] 



S = 



A = 



Ccj) 



(14) 



which coincide with (|l3|), (p"2|), (]?]) after ignoring sub- 
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and t p for numerical solution of 



TABLE I: Parameters 

Eqs. ( pi] ) - ( |l3] ) with C\ = C2 = \, and computed values of 
S, crossover scale, k», and normalized rates of dissipation: by 
the iteration procedure, /200, /300 and numerically, f^. 






1.0 


0.75 0.5 0.25 




1.0 0.5 0.2 0.1 0.02 


0.5 


8 


0.21 0.10 0.04 0.02 0.004 
0.42 0.20 0.08 0.04 0.008 
3.4 9.5 39 117 1530 


0.12 0.15 0.19 
0.21 0.22 0.24 
9.2 8.8 8.4 


j3,oo 
foe 


0.704 0.621 0.522 0.455 0.331 
0.721 0.639 0.537 0.469 0.341 
0.723 0.642 0.541 0.474 0.352 


0.686 0.766 0.868 
0.699 0.774 0.871 
0.701 0.775 0.871 
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FIG. 1: Log-Log plots of solutions of Eqs. (0) - @ for 
parameters in Tab. |j. Solid lines: "exact" numerical solutions 
f(kL) = e{k)/e L , dashed lines: f 2 {kL) = e 2 (k)/e L , Eq. @. 



script "„" . First two iterations may be done explicitly 

El 

Peff,l(fc) = 1 



ti(k) 
/a(«) 



ffl(K) = 2lW =K (^A)/3. 



7l 



p f 1 + t p 7 i _k( 2 - a )/ 3 

ea(*) _ r<5 + K ( A - 2 )/ 3 i3A/(2-A) 



(15) 
(16) 



S + l 

g 2 ( K ) = =« 2/3 / 2 1/3 (^K 1/3 W, 

r2(t) = p c s,2{k)/pf = l + 0/[l + Tp7i 52 («)] • 
The third iteration requires simple numerical integration: 

e 3 (fc) =e i exp| - A J dx/x [1 + ^(x)] } . (17) 

Figure [l] displays "exact" numerical solutions of 
e(k)/e L (for sets of </> and r p listed in Tab. |) in com- 
parison with corresponding plots of Eq. (Jig) for fzfa). 



In the first decade of the inertial interval (kL < 10) dis- 
crepancies between the numerics and the second itera- 
tion are quite small. They gradually increase toward 
large k, remaining smaller than (1 — 2)%. Values of 

/„, = /(oo), / 2 ,oo = [5/(1 + 5)} 3A/(2 - A) and / 3)00 fol- 
lowing from Eq. ( |i~7| ) are given in Tab. |. It is clear that 
almost always one may use simple analytical solution of 
the second iterations, dH). The results @ of the third 
iteration may be used to control accuracy. 
3. Turbulence modification by particles. In the 
particle-free turbulence the rate of energy input e L is 
equal to the energy flux in the inertial interval e(k) 
and to the rate of energy dissipation e^. In turbu- 
lent suspensions due to fluid-particle friction e(k) is no 
longer constant and decreases toward large k. Therefore 
e L > e(k) > Soo. Eq. ( p"5| ) and Fig. |l]a show that in the 
region k <C fc* [ T pd(k*) = 1] the flux decays toward small 
scale: e(k) « e L (kL)~ A with A, Eq. (B), independent of 
r p , while for k 3> fc* the flux approaches plato, similar to 
the particle-free case. At first glance these two facts are 
unexpected: there is an essential suppression of the large 
scale eddies in spite of the fact that particles are almost 
swept by them and therefore the fluid-particle damping 
7 P (k) is small. On the contrary, the small scale eddies are 
almost not effected by particles which are not involved in 
their motion and thus 7 P (fc) reaches the maximum value 
4>/t p . To explain consider the dimensionless rate of the 
flux decrement, [— dhie(k) /dink], which is ex 7 P (fc). The 
only available frequency to make 7 P (fc) dimensionless is 
7(fc). Therefore -dlne(k)/dlak ~ T(k) = 7 P (fc)/7(fc), 
in agreement with Eq. (|l3|). As one sees from Eq. (0) 



r(k) = <j>/[i + <p + T pl (k)] 



(18) 



and for r p 7(fc) <C 1 the ratio T(k) becomes t p indepen- 
dent constant. This explains both facts: a - why the sup- 
pression of small k eddies is t p independent and b - why 
this suppression is scale invariant. Note that the weak 
sensitivity of turbulent spectra to t p (fact a) was previ- 
ously observed in numerous simulations of turbulence in 
dilute suspensions but, to the best of our knowledge, was 
not well understood, see, e.g. ||. Equation ( |l8| ) shows 
that r(fc) — > for r P 7(A;) 3> 1. Therefore in this region 
of k particles cannot modify the turbulence and indeed 
e(k) must approach a constant value, 600 . 

For brief comparison of the prediction b with direct nu- 
merical simulation by Boivin, Simonin and Squires fi) we 
reploted in Fig. || their Fig. 5b for kinetic energy spectra 
E(k) of suspensions in Log-Log coordinates (lower lines). 
Our first iteration ( |l5| ) predicts E(k) oc fc-5/3-2A/3 w ^ 
A = C<p/(l + 4>). Therefore we defined "compensated" 
spectra as E c {k,cj>) = E{k)(kL) 5 / 3+2A / 3 and plotted 
them in Fig. | (with C = 3.8). The region 0.4 <Log 
k < 1 where for = 0, E c (k) = E(k)(kL) 5 / 3 , upper solid 
line, is approximately constant may be considered as in- 
ertial interval. As we expected, in this interval all lines 
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FIG. 2: Log-Logplots of turbulent kinetic energy spectrum 
E(k) taken from M for cf> = 0, 0.2, 0.5 and 1 (lower lines) and 
"compensated" spectra E c (k) for the same <j) (upper lines). 
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FIG. 3: Upper panel: Log-log plots of functions J?2(fc) = 
[/2(K)r 2 (l)/r 2 (K)] 2/3 , Eqs. @, (!§), for = 1 and values of 
r p denoting corresponding lines. Lower panel: Plots of 7?2,oo, 
Eq. (|20[), vs. r p . Values of label corresponding lines. 



for different <f> are about to collapse. Some scattering of 
the lines is related with fc-dependence of p e s{k) and fi- 
nite value of r p which is neglected in the first iteration. 
More detailed comparison of our second and next order 
iterations will be done elsewhere [fjj. 

Consider now a possible enhancement of the density of 
kinetic energy of the carrier fluid E(k). According to Eqs. 

(§), <§ and @ E(k) = C m [e(k)/ PcS (k)] 2/3 k- 5 / 3 . In- 
troduce the dimensionless ratio 



R{k) 



E(k)/E(l/L) 
E (k)/E (l/L) 



e(k)p L 



-i 2/3 



(19) 



Here E (k) = Ci£ 2/3 p f 1/3 fc- 5/3 is the density of turbu- 
lent kinetic energy in the particle-free case. This ratio 



is larger (smaller) than unity in the case of enhancement 
(suppression) of the turbulent energy by particles. To 
understand this behavior consider three distinct regions 
of scales defined by the crossover scale £;* ( see Table |) : 
a. Region of small scales, L^ 1 < k < k*, where e(fc) is 
decreasing function of fc. Function p c s(k) in the denom- 
inator of Eq. ( |l9| ) is almost constant pf(l + </>). 

b. Region of transient scales, k ~ fc*, where e{k) still 
decreases, and so does p e s(k) gradually reducing to pf. 

c. Region of large scales, k > > rj~ l , where e{k) 
approaches plato, while p e s(k) is again constant, pf. It 
is clear that behavior of R(k) will depend on k*L. For 
k*L 3> 1 (small enough r p ) the ratio R(k) has enough 
room in the region a to strongly decay [as (kL)~ 2A ^ 3 ] 
due to decay of e(k). In the region b it may increase 
[due to decrease of p c g(k)] but not more than by factor 



of (1 



)2/3 



Therefore for small enough t p the ratio 



R(k) < 1 everywhere, see, e.g., in Fig. ||a plots for t v — 
0.05, ( KL « 350) and t p = 0.1 (k*L = 117). 

There is an essential enhancement of the kinetic energy 
for larger r p (smaller k^L), when the small scale region 
a is not pronounced (plots for t p = 0.5 (k^L = 9.5) 
and r p = 0.75 (k^L — 5.2) in Fig. ||a). In this case the 
growth of R(k) in the transient region b is stronger than 
the decay in the region a. 

To find values of parameters for which the enhance- 
ment is possible (at least for k — > oo) consider the maxi- 
mum value of i?oo = R{po). In the second iteration 



-R2,oo — 



1 



2A/(2-A) 



1 



5<p 



1+5 



2/3 



(20) 



Fig. ||b displays plots i?2,oo vs. r p for different <f>. The 
enhancement is possible if 8 = t p 7^/(1 + <p) exceeds a 
critical value <5 cr which is independent of <f>. Parameter 
cj) governs the value of the enhancement. The maximal 
possible enhancement (for large kL and S) is (1 + (j)) 2 ^ 3 - 
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